Gap soliton dynamics in an optical lattice as a parametrically driven pendulum 
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A long wavelength optical lattice is generated in a two-level medium by low-frequency contraprop- 
agating beams. Then a short wave length gap soliton generated by evanescent boundary instability 
(supratransmission) undergoes a dynamics shown to obey the Newton equation of the parametrically 
driven pendulum, hence presenting extremely rich, possibly chaotic, dynamical behavior. The theory 
is sustained by numerical simulations and provides an efficient tool to study soliton trajectories. 
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Introduction. Nonlinear physics has revealed that 
quite different complex systems may actually share the 
same model equations with common simple and univer- 
sal physical properties [1]. One celebrated example is 
the Fermi-Pasta-Ulam chain of anharmonic oscillators 
[2] , which may serve as a laboratory to check soliton the- 
ory, statistical physics and even dynamical processes in 
DNA molecules. Another famous example is the Joseph- 
son junction, mathematically analogue to a pendulum, 
where the biased voltage and the internal resistance play 
the role of forcing and damping [3[. 

We consider here another well established Maxwell- 
Bloch (MB) model Q which has an universal character 
related to various nonlinear processes [5|, Q in nonlin- 
ear optics, and discover a new example of a reduction 
of a complex many body dynamics of MB system to the 
driven-damped pendulum motion. This is done in the 
context of gap soliton motion in a two level medium sub- 
jected to some low frequency stationary boundary driv- 
ing. We show that the soliton trajectory in this effec- 
tive optical lattice (periodic in space and time) obeys 
the equation of a parametric pendulum which results in 
a rich dynamical behavior, from periodic to chaotic, de- 
ending on the fundamental parameters of the problem 
1[. Fig H] displays three instances of the gap soliton 
propagation in a two level medium prepared as an op- 
tical lattice. The trajectories compare well to the time 
dependence of the angle of a parametric pendulum. 

In a two-level system (TLS) of transition frequency fio, 
the governing equation is the MB model [4|], considered 
here in the isotropic case for a linearly polarized electro- 
magnetic field propagating in direction z. The time is 
scaled to the inverse transition frequency Qq 1 , the space 
z to the length VLqc/t] with optical index n = ^/e/i , the 
energy to the average Wo = iVoftfio/2, the electric field 
to y/Wo/e and the polarization to y/eWo- Here No is the 
density of active dipole. The resulting dimensionless MB 
system then reads 

(E + P) tt = E zz , P u + P= -aNE, N t = EP t . (1) 

The electric field is E = (E(z, £), 0, 0) and the polariza- 
tion source P = (P(z, £), 0, 0). The quantity N(z,i) is 



the normalized inversion of population density that is 
assumed to be —1 in the fundamental state (no applied 
field). The coupling strength is the dimensionless funda- 
mental constant a completely characterized by the gap 
opening between the transition frequency (value 1 in re- 
duced units) and the plasma frequency c^o- This gap 
results from the linear (TV = —1) dispersion relation of 
the MB system 



uo 2 (uo 2 -uo 2 ) = k 2 (co 2 - 1), 



uj 2 



(2) 



obtained for a carrier ex.-p[i(cot — kz)]. 

Our purpose is to create an optical lattice with two 
contra-propagating beams of frequency Q < 1, which 
then will interact with a wavepacket of central frequency 
w < uoq (inside the gap, close to the upper gap edge). 
The contra-propagating beams are expected to create a 
stationary wave in the variables E, P and TV that will 
then interact with the gap soliton through mediation of 
the two fundamental coupling terms NE and EP t in MB 
equations. To achieve this study we shall derive from 
(PQ) a limit model by the reductive perturbative expan- 
sion method where essential phase effects are carefully 
taken into account [9|. We shall obtain the following dy- 
namical equation for the soliton motion, in some newly 
normalized time t and position q(t): 
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(3) 



which is a parametrically driven pendulum. Here £ is 
the amplitude of the stationary standing wave and 7 is a 
small phenomenological damping parameter accounting 
for soliton energy losses trough the optical lattice. The 
comparison of this effective dynamics with the numerical 
simulations on the full MB model is presented in fig{H 

Theory. We consider an electric field that carries two 
fundamental frequencies: one is close to the band gap 
edge cjo, the other is close to zero. Then one deals with a 
two-wave nonresonant interaction process whose weakly 
nonlinear limit is sought by assuming the formal series 
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FIG. 1: (Color online) Graphs (a), (b) and (c) display results 
of numerical simulation of (Jj) under boundary conditions Offi) 
for the three incident wave amplitudes (24\) generating opti- 
cal lattice with slightly different depths. Graph (d) shows the 
pendulum angle evolution in model {2^ corresponding to those 
three regimes. For sake of presentation, the gap soliton tra- 
jectories have been underlined with dashed lines. 



expansion [9[ 



p—1 £=—p 

^=s[z-^ 2 ,r)}, ti=s 2 z, r = e% (4) 

where F stands for any of the three fields E(z, t), P(z, t) 
or n(z, t) — l-\-N(z, t). Note that now n(z, t) denotes the 
normalized population of the excited state. As all com- 



ponents are real valued, we have F^ = F^ J (over bar 
stands for complex conjugate). The slow space variables 
£1 and £2 are associated respectively with the character- 
istic wavelengths of the gap soliton and of the standing 
wave grating. The above representation actually means 
to replace the differential operators as follows 
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dz 
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dt 
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(5) 



when applied to a given harmonic £ of the asymptotic 
series inserted in (pQ). This provides the following first 



order structure of the field 



P = -e[E[ 1 \^r)e 



c.c.l 



n=0(e 2 ), 



(6) 



which accounts for the linear limit and deserves some 
comments. The set of unknowns is: the slowly varying 
envelope £^(£i,r) of the short wavelength wavepacket 
(e.g. gap soliton), the very slowly varying profile 
4^(6j t ) °f the long wavelength applied grating, the 
phase #(£2? T ) which accounts for variations due to wave- 
coupling and finally ^(£2, t) which will be related to the 
position of the wavepacket by means of (|4j). Then we 
note that the expression of the polarization in terms of 
the electric field actually does follow the usual linear laws 
P — —E for the plasma wave (here given by E^ e luJot , 
the mode at uo — luq) and P = aE for the electrostatic 
component (here given by E^\ the m °de at u = 0). 
The dependence of these two fields in the slow variables, 
together with phase variations, will then carry the non- 
linear electrodynamics. 

Next we seek the leading order (e 2 ) harmonics of 
n(z,t). According to the last equation of (pQ), using ([6]) 
and collecting harmonics i — 1 and i = 2, we get after 
time-integration 



,(2) 



>1 



0)\ 



(7) 
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The last term to compute is which cannot be calcu- 
lated at order e 2 since there, the third equation of (JT|) is 
trivial. Thus we move one step further (actually to e 4 to 
catch the r-derivative), which provides 
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(8) 



Then writing the second equation of (pQ) at order e 2 and 
using 0, we obtain for i = 1 



a(P (3) +£ ;(3)) + 2^ ^-2^ 1 )g 



(9) 



4 + Q lr7(l)l 2 i7(l) 



E™\'E™ + a(E^) 2 E^ - cmg'E™ = 0, 

which is replaced in dSj). The result can be integrated to 
furnish the sought expression 
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Next we consider first equation of (pQ) at order e 3 and 
collect again the £ = 1 harmonics to get 



(3) 



(1) 



0. 



After some simple algebra we can eliminate E^ and P| 
from equations ([9]) and (fTTj). The resulting equation ap- 
pears with terms that depend solely on a single variable, 
e[^ (£i ? t) on the one side, 0(^2 •> t) and (£2 , t) on the 
other side. These terms thus decouple to eventually give 



(11) 



,(3) 



dE[ 
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(12) 
(13) 



It appears therefore that E^\^i,r) obeys a nonlinear 
Schrodinger equation where the effect of the applied grat- 
ing Eq 1 ^ lies in the definition of the variable £1 in (|4|). 
Thus the drift 4>{^2 jt) remains to be evaluated, which is 
done with the last equation of (pQ) at order e 4 and i = 1 



dr 



a 86 
^0 d& 



(14) 



(2) 

Here we have assumed E\ = which is allowed by the 
structure of ([1]) where nonlinearity comes into play at or- 
der e 3 . Indeed one actually obtains a linear homogeneous 



equation for e[ 2 ^ whose unique solution is E\ A) = as 
soon as the initial-boundary value problem concerns only 

To close the system ([T2]) - ([H]) we derive the equations 
for the low frequency oscillations obtained from the first 
equation of (pp) at £ = 0, namely 



(2) 



,2 d 2 E 



(i) 



8 2 E ( 



(i) 



dr 2 



d& 2 



= 0. 



(15) 



The system ([T2]) - ([T5]) is the basic set of equations that de- 
scribes within the Maxwell- B loch model, the interaction 
of a gap soliton with a long wavelength standing wave. 
Note that all coefficients are completely determined from 
the unique parameter a, the fundamental coupling con- 



stant of eq.(pp), as by definition 1 
eliminate the phase between ([I 

d 2 <\> 
dr 2 



- 1 + a. It is useful to 
and (p3j) and obtain 



a 2 (a-2) d 
4(1 + a) 2 d& 



(4 



(16) 



which constitutes with ([T5]) a closed system, indepen- 
dently of E^\ 

Application. We apply now the above machinery to 
the case when the fundamental field component is 
a gap soliton, exact solution to ([T2]h of given velocity 
in the frame (£i,t). The nature of equation ([T2]) allows 



us to start with such an explicit solution and then to 
study its dynamics by looking at the frame drift ^(£2, t) 
which thus defines the variations of the soliton about its 
free motion. To that end it is more convenient to rewrite 
system ([T5]) ([T6 ]) in the physical (dimensionless) variables, 
namely 



2^0 

d 2 



d^Eo 

dt 2 dz 2 
a 2 (a-2) d 



0, 



dt 2 4(1 + a) 2 dz 



(17) 



where E (z,t) = eE^^r). 



Now the soliton is lo- 
say z s (t), is defined by 



calized in £1 and its position 
£1 = const, namely by z s — <j)(z Sl t) = const. Computing 
the velocity dz s /dt and the acceleration d 2 z s /dt 2 (total 
derivatives) we remember that <j) is slowly varying in z 
and t and thus keep only the dominant orders to eventu- 
ally obtain 



d 2 Z 



dt 2 



dt 2 



(18) 



We may now select a particular solution to the linear 
wave equation ([T5]) that would result from application 
to the medium of two contrapropagating monochromatic 
beams. The effective boundary conditions that represent 
such a situation will be described later, they result in the 
generation of a standing wave at frequency namely 
(assuming a system length corresponding to a mode at 
that frequency) 



E = £ sin(ftt) sm(Kz), K 2 = (1 + a)n 2 



(19) 



Note that the above dispersion law is indeed the behavior 
at small uo of the general dispersion law ([2]) . The resulting 
equation for the soliton acceleration ([15)1 then reads 



d 2 z s 
dt 2 



KS 



a 2 (a - 2) 
4(1 + a) 2 



sin 2 (f}t) sm(2Kz s 



(20) 



This is the parametric driven pendulum equation that 
can be written as ((3|) by rescaling time (fit — > t) and 
position (2Kz s = q), and by adding a phenomenological 
damping to account for soliton radiation. Such an equa- 
tion is solved with the data of initial soliton position zo 
and velocity vq. Our purpose now is to compare the so- 
lution of Maxwell-Bloch (pQ) with convenient boundary 
data to the prediction of soliton trajectory given by ((3|). 

Numerical simulations. To proceed with numerical 
simulations of the MB equations (pQ) we derive the 
boundary- value data following [ll|, that represent inci- 
dent waves entering the system at z = and z = L 
and which are expected to produce the standing low fre- 
quency wave ([T9|) . The vacuum outside [0,L] is assumed 
to obey (pQ) with a — whose solution reads 

z < : E vac = Icos[fl(t - z)] + Rcos[ft(t + z)}. (21) 



4 



The amplitude / of the incident wave from the left is the 
control parameter, the amplitude R of the reflected wave 
has thus to be eliminated. The continuity conditions at 
z = with the electric field E(z, t) inside the medium 
can be written 

(I + R)cos(nt) = (E) z=0 , 
n(I-R)sm(nt) = {d z E) z=Q , 

which can be combined to eliminate the unknown re- 
flected amplitude R to give 

(d z E - d t E) z=Q = 20Jsin(Ot), 
{d z E + d t E) z=L = -2QQsin(ft£). 

The above second relation is obtained at t z = L and 
Q is the chosen amplitude of the incident wave from the 
right. Such unusual boundary data are actually numeri- 
cally implemented in a finite difference scheme as (define 
L~ =L- dz) 

E(0, t) =E(dz,t)- dz [d t E(dz, t) + 20/ sin(ftt)] , 
E(L,t) =E(L~,t) - dz [dtE(L~,t) + 2fiQsin(fit)], 

which are inserted in the differential equation (pQ). We 
shall use equal amplitude out of phase driving, namely 
I = —Q. This generates the standing wave solution (Qi 
whose amplitude £ is then defined from / by 



£ = 21 [(1 + a) cos 2 (KL/2) + sin 2 (KL/2) 



-1/2 



. (22) 



Then the gap soliton is generated by driving the 
boundary z = with an incident pulse at frequency uo 
[Tot and finally the full set of boundary conditions that 
produce the plots (a), (b) and (c) in fig H] read 

(d z E-d t E) = 2I0sin(0t) + 

v Jz=0 v ; cosh[(t- 100)/12]' 

(d z E + d t E) 2=L = 2m sin(Qt) , (23) 

for a system length L = 254 and initial state at rest: 
E(z,0) = 0, P(z,0) = 0, N(z,0) = -1. We choose 
a = 1, Q = tt/200 and u = 1.5. As defined by ( ED the 
generated standing wave has wavenumber K = Q^l + a 
and its amplitude is defined by ([22]) . The only parameter 
we vary is the amplitude / of the contrapropagating low 
frequency beams. In particular the graphs (a), (b) and 
(c) of figtU correspond to the following values 



2I a = 1.832, 2I b = 1.843, 



2L 



1. 



(24) 



Why such a small change of the control parameter causes 
completely different trajectories of the soliton is under- 
stood using the pendulum description (|3]) where in view 
of (|22|) the amplitude A is related to the parameters of 
MB (with a = 1) by 



We now solve (|3j) with initial conditions z s (0) = 0.8 and 
= 0.94 that have been measured as the common gap 
soliton initial position and velocity in the three succes- 
sive simulations (a), (b) and (c). The phenomenological 
damping constant is 7 = 0.01. Then the pendulum angle 
evolution for three different values of A resulting from 
([25]) and ([24]) are displayed in the graph (d) of figffl 

Conclusion and comments. It is remarkable that such 
different models as the Maxwell-Bloch system (pQ) and 
the driven parametric pendulum ([3]) concur to describe 
the dynamics of a gap soliton interacting with a low- 
frequency standing wave. This demonstrates in particu- 
lar that the gap soliton trajectory may well be chaotic, 
at least initially (it may stabilize due to energy losses by 
radiationor the internal damping of MB model). This is 
a direct result of the time dependence of the grating in- 
duced by the stationary field. In particular the gap soli- 
ton dynamics in a permanent grating (as e.g. a Bragg 
medium) is regular. 

It is worth noting that the effect of damping in initial 
MB model (pQ), originating from finite dephasing times, 
causes eventual transition from chaotic to self trapped 
regime, as expected from a driven damped pendulum, 
and which we have checked on numerical simulations. 
This occurs for time scales of the order of the dephasing 
times, for which the soliton gets trapped into one of the 
nodes of the standing wave grating. 

One practical advantage of the parametric pendulum 
description is the prediction of the switching from regular 
to chaotic (erratic) motion of the gap soliton in terms of 
the parameters of the driving field. However we expect 
that the various scenario of chaotic versus periodic mo- 
tion could be more complex than the driven pendulum 
case and a more comprehensive analysis would first re- 
quire much longer computation times, which is a subject 
of further studies. Beyond such a technical result, our 
analysis has revealed a new type of soliton motion in a 
medium with dynamically generated optical grating. 

We expect to extend this analysis of a one dimensional 
time-dependent situation to the study of spatial soliton 
trajectories in two dimensional smooth optical lattices. 
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A = I 2 \2cos 2 (KL/2) + sin 2 (KL/2) 



(25) 
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